﻿#region Copyright information
// 
// Copyright © 2005-2013 Yongkee Cho. All rights reserved.
// 
// This code is a part of the Biological Object Library and governed under the terms of the
// GNU Lesser General  Public License (LGPL) version 2.1 which accompanies this distribution.
// For more information on the LGPL, please visit http://bol.codeplex.com/license.
// 
// - Filename: SelfOrganizingMap.cs
// - Author: Yongkee Cho
// - Email: yongkeecho@gmail.com
// - Date Created: 2012-12-08 8:26 AM
// - Last Modified: 2013-01-25 3:59 PM
// 
#endregion
using System;
using System.Collections.Generic;
using System.Linq;
using BOL.Linq.Descriptive;

namespace BOL.Algorithms.Clustering
{
    /// <summary>
    /// Implements a self-organizing map clusterer.
    /// <see>http://en.wikipedia.org/wiki/Self-organizing_map</see>
    /// <see>http://www.cis.hut.fi/somtoolbox/documentation/somalg.shtml</see>
    /// </summary>
    public class SelfOrganizingMap : ClustererBase<double>, IMapClusterer, IIterative
    {
        #region Private variables

        private readonly Func<double, double> _neighborhoodFunction; 
        private readonly double _initialLearningRate;
        private readonly double _finalLearningRate;
        private int _gridCount;

        #endregion

        #region Public properties

        public double[][] Map { get; private set; }

        public int GridSize { get; private set; }

        public int NumberOfIterations { get; set; }

        #endregion

        #region Constructor

        public SelfOrganizingMap(int gridSize, Func<IEnumerable<double>, IEnumerable<double>, double> distanceFunction, Func<double, double> neighborhoodFunction, int numberOfIterations = 100, double initialLearningRate = 0.9, double finalLearningRate = 0.1)
            : base(distanceFunction)
        {
            GridSize = gridSize;
            _gridCount = gridSize * gridSize;
            
            _neighborhoodFunction = neighborhoodFunction;
            NumberOfIterations = numberOfIterations;
            _initialLearningRate = initialLearningRate;
            _finalLearningRate = finalLearningRate;
        }

        public SelfOrganizingMap(Func<IEnumerable<double>, IEnumerable<double>, double> distanceFunction, Func<double, double> neighborhoodFunction, int numberOfIterations = 100, double initialLearningRate = 0.9, double finalLearningRate = 0.1)
            : this(0, distanceFunction, neighborhoodFunction, numberOfIterations, initialLearningRate, finalLearningRate)
        {
        }

        #endregion

        #region Private methods

        private int FindBestMatchingUnitIndex(IEnumerable<double> v)
        {
            return Map.Select(u => DistanceFunction(u, v)).ArgMin();
        }

        #endregion

        #region Public methods

        public override IEnumerable<IGrouping<int, TSource>> Cluster<TSource>(IEnumerable<TSource> sources, Func<TSource, double[]> pointSelector)
        {
            if (sources == null)
                throw new ArgumentNullException("sources");
            if (pointSelector == null)
                throw new ArgumentNullException("pointSelector");

            var list = sources.ToList();
            var points = list.Select(pointSelector).ToList();
            var numberOfPoints = points.Count;
            if (numberOfPoints == 0)
                throw new ArgumentOutOfRangeException("sources");
            var length = points[0].Length;
            
            // Determines grid size if undetermined or inappropriate
            if (GridSize < 1)
            {
                GridSize = (int)Math.Sqrt(2 * numberOfPoints);
                _gridCount = GridSize * GridSize;
            }

            // randomizes map nodes' weight by sampling within the data range
            var r = new Random();
            var range = Enumerable.Range(0, length).Select(i => points.Range(x => x[i])).ToArray();
            Map = new double[_gridCount][];
            Map = Map.Select(x => new double[length]).ToArray();
            for (var j = 0; j < length; j++)
                for (var i = 0; i < _gridCount; i++)
                    Map[i][j] = r.NextDouble() * range[j];

            // iterates till the max
            var iter = 0;
            do
            {
                // learning rate decays polynomial
                var learningRate = _initialLearningRate * Math.Pow(_initialLearningRate / _finalLearningRate, (double)iter / NumberOfIterations);

                // for each data point,
                foreach (var point in points)
                {
                    // finds the best matching unit
                    var bmu = FindBestMatchingUnitIndex(point);

                    // for each map,
                    for (var i = 0; i < _gridCount; i++)
                    {
                        // estimates circular distnace
                        // TODO: implement rectangular or hexagonal distance here
                        var distance = Math.Sqrt(Math.Pow(i / GridSize - bmu / GridSize, 2) + Math.Pow(i % GridSize - bmu % GridSize, 2));

                        // updates weights 
                        Map[i] = Map[i].Zip(point, (m, p) => m + _neighborhoodFunction(distance) * learningRate * (p - m)).ToArray();
                    }
                }
            } while (iter++ < NumberOfIterations);

            // maps data point into the cloest map point
            return list.GroupBy(x => FindBestMatchingUnitIndex(pointSelector(x)));
        }

        #endregion
    }
}
